// Couper une cyclide par des plans correspondant aux cercles de Villarceau
// On essaie de rendre le tout un peu plus analytique

import ij.IJ;
import ij.ImagePlus;
import ij.gui.NewImage;
import ij.io.FileSaver;
import ij.plugin.PlugIn;
import ij.process.ColorProcessor;
import ij.process.ImageProcessor;


import java.awt.Color;


import java.awt.Font;

import java.awt.Polygon;
import java.io.BufferedReader;
import java.io.DataInputStream;
import java.io.FileInputStream;
import java.io.FileReader;
import java.io.PrintStream;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.Locale;
import java.util.Random; 
import java.util.Scanner;



import pluginUtils.Circle;
import pluginUtils.Circle3D;
import pluginUtils.CircleArrangement;
import pluginUtils.Complexe; 
import pluginUtils.DoyleCircle;
import pluginUtils.Droite;
import pluginUtils.MobiusTransform;
import pluginUtils.CircleInversion;
import pluginUtils.PackCircle;
import pluginUtils.Plane;
import pluginUtils.Pos3D;
import pluginUtils.Sphere;
import pluginUtils.SphereInversion;
import sliceformUtils.CoupleDAngles;
import sliceformUtils.DroiteSpace3D;
import sliceformUtils.Ligne;
import sliceformUtils.Mesh;
import sliceformUtils.STLReader;
import sliceformUtils.Saucisson;
import sliceformUtils.Segment;
import sliceformUtils.Triangle;


public class SliceFormCyclideV3 extends ImagePlus implements PlugIn{
	{Locale.setDefault(Locale.US); }
	
	private ImageProcessor primeIP; 
	private int X=2000,Y=2000; // image size
	
	private DataInputStream input;
	private double xShift=0; 
	private double zShift=0; 
	
	private static int nbCercles=6; 
	
	private static Circle3D lesCercles[]=new Circle3D[nbCercles];
	
	
	private double a=5,mu=3,c=1;
	
	Pos3D pointCommun=new Pos3D(a*c/mu,0,0); 
	
	
	
	private DroiteSpace3D intersections[][];
	
	

	private Random generator=new Random(); // One always needs some randomness
	
	// Scaling  the image
	private double x0=0,x1=1,y0=0,y1=1;
	private double acoef,ccoef,vcoef,wcoef; 
	
	public void setCoefs(){
		acoef=X/(x1-x0); 
		vcoef=Y/(y1-y0); 
		ccoef=-acoef*x0; 
		wcoef=-vcoef*y0; 
	}
	 
	
	
	public void setBounds(double xl,double xu,double yl,double yu){
		x0=xl; 
		x1=xu;
		y0=yl;
		y1=yu; 
	}
	
	private void drawCircle(Circle3D c,ImageProcessor i1){
		double xCenter=c.getCenter().getY(); 
		double yCenter=c.getCenter().getZ(); 
		//if(c.getRadius()<0.001) return ; 
		
		int intX=(int)((xCenter-xShift)*acoef+ccoef);
		int intY=(int)((yCenter-zShift)*vcoef+wcoef);
		int intRx=(int)(c.getRayon()*acoef);
		int intRy=(int)(c.getRayon()*vcoef);
		
		
		i1.drawOval(intX-intRx,intY-intRy,intRx*2,intRy*2);
	}
	
	private void drawPoint(double x,double y, ImageProcessor ip){
		int intX=(int)(x*acoef+ccoef);
		int intY=(int)(y*vcoef+wcoef);
		ip.fillOval(intX-10,intY-10,20,20);
	}
	
	private void drawLine(double xd,double yd,double xf,double yf){
		int XD=(int)(xd*acoef+ccoef); 
		int YD=(int)(yd*vcoef+wcoef); 
		int XF=(int)(xf*acoef+ccoef); 
		int YF=(int)(yf*vcoef+wcoef);
		primeIP.drawLine(XD,YD,XF,YF); 
	}
	
	private void drawLine(double xd,double yd,double xf,double yf,ImageProcessor ip){
		int XD=(int)(xd*acoef+ccoef); 
		int YD=(int)(yd*vcoef+wcoef); 
		int XF=(int)(xf*acoef+ccoef); 
		int YF=(int)(yf*vcoef+wcoef);
		ip.drawLine(XD,YD,XF,YF); 
	}
	
	// Couper un Mesh par un plan
	public Ligne tranche(Mesh leMesh,Plane lePlan){

		Ligne contour=new Ligne(lePlan);
		Iterator<Triangle> it=leMesh.iterator(); 
		while(it.hasNext()){
			Triangle tr=it.next(); 
			if (tr.isCut(lePlan)) {  contour.add(tr.coupePar(lePlan));}		
		}
		
		//System.out.println("\n Contour :\n"+contour.info()); 
		
		
		return contour;
	}
	
	
	// b=true : on affiche les P1, b=false, on affiche les P2 (utile pour tracer les lignes de coupures)
	public void print(Ligne olida,Circle3D circles[],double wide,double marge,String name,int numero,boolean couleur){
	    name="dupin/dir3/YYY"+name;
		double bmin=olida.getMinY();
		double cmin=olida.getMinZ();
		double bmax=olida.getMaxY();
		double cmax=olida.getMaxZ();
		double centerY=(bmax+bmin)/2; 
		double centerZ=(cmax+cmin)/2; 
		System.out.println("Bornes mins : "+bmin+" "+cmin+"\n"); 
		setBounds(centerY-marge-wide/2,centerY+marge+wide/2,centerZ-marge-wide/2,centerZ+marge+wide/2); 
		setCoefs();
		
		
		CoupleDAngles larot=rotations(circles[numero].getPlan());
		
		ImagePlus gamma=NewImage.createRGBImage (name+numero,X,Y, 1, NewImage.FILL_WHITE);	
		ImageProcessor i1=gamma.getProcessor();
		FileSaver fs=new FileSaver(gamma); 
		i1.setColor(Color.RED);
		i1.setLineWidth(5);
		drawLine(centerY-wide/2,centerZ-wide/2,centerY+wide/2,centerZ-wide/2,i1); 
		drawLine(centerY-wide/2,centerZ-wide/2,centerY-wide/2,centerZ+wide/2,i1); 
		drawLine(centerY-wide/2,centerZ+wide/2,centerY+wide/2,centerZ+wide/2,i1); 
		drawLine(centerY+wide/2,centerZ-wide/2,centerY+wide/2,centerZ+wide/2,i1); 
		
		Pos3D pointCourant=pointCommun.rotate(larot);
		drawPoint(pointCourant.getY(),pointCourant.getZ(),i1);
		
		i1.setLineWidth(4);
		if(couleur)
			i1.setColor(Color.BLACK);
		else
			i1.setColor(Color.BLUE);	
		
		Iterator<Segment> it=olida.iterator();
		while(it.hasNext()){
			
				Segment seg=it.next();
		
				
				Pos3D s1=seg.getOne(); 
				Pos3D s2=seg.getTwo(); 
				drawLine(s1.getY(),s1.getZ(),s2.getY(),s2.getZ(),i1); 
			}
		
		//toutes les lignes
		for(int j=0;j<nbCercles;j++){
			if(j!=numero){
				DroiteSpace3D courante=intersections[j][numero]; 
				System.out.println("Droite avant rotation : "+courante); 
				courante=courante.rotate(larot); 
				System.out.println("Droite apres rotation: "+courante); 
				
				ArrayList<Pos3D> lesBords=intersectCadre(centerY-wide/2,centerY+wide/2,centerZ-wide/2,centerZ+wide/2, courante);
				if(lesBords.size()==2){ // tracer la ligne
					Pos3D s1=lesBords.get(0); 
					Pos3D s2=lesBords.get(1);
					//System.out.println("--->"+s1); 
					drawLine(s1.getY(),s1.getZ(),s2.getY(),s2.getZ(),i1); 
				}// if lesBords
				
				
				
			}
			
		}
			
			// centre du cercle DEBUG
			drawPoint(lesCercles[numero].getCenter().getY(),lesCercles[numero].getCenter().getZ(),i1);
			i1.setLineWidth(4);
			i1.setColor(Color.RED) ;
			//drawCircle(circles[numero],i1);
			drawCircle(lesCercles[numero],i1);
			gamma.show(); 
			fs.saveAsPng("/tmp/"+name+numero+".png"); 
			//gamma.close(); 
			
			
		}

		
	
	
	/** Calcule les angles de rotation amenant un plan parallele a un plan perpendiculaire a ox **/
	
	public CoupleDAngles rotations(Plane p){
		CoupleDAngles resu=new CoupleDAngles(); 
		if(p.getA()==0){
			if(p.getB()==0){
				resu.setAlpha(0); 
				if(p.getC()<0) resu.setBeta(Math.PI/2); else resu.setBeta(-Math.PI/2); 
				//System.out.println("resu "+resu); 
				return resu;
			}
			else // A=0 et B!=0
			 {
			 if(p.getB()<0) resu.setAlpha(Math.PI/2); else resu.setAlpha(-Math.PI/2);
			 }
			 }
	    else // A!=0; 
	    	resu.setAlpha(-Math.atan2(p.getB(),p.getA())); 
			
		// On s'occupe de beta
		double leng=Math.sqrt(p.getA()*p.getA()+p.getB()*p.getB()); 
		resu.setBeta(-Math.atan2(p.getC(),leng));
		
		//System.out.println("resu "+resu); 
		return resu; 
		}
		
	public DroiteSpace3D inter(Plane p1,Plane p2){
		if((p1.getA()==0)&&(p2.getA()!=0)){
			Plane temp=p1; p1=p2; p2=temp;
		}
		if((p1.getA()==0)&&(p2.getA()==0)){ // droite de direction (1,0,0) (ou bien plan paralleles, pas de solution)
			if((p1.getB()==0)&&(p2.getB()==0)) return null; // plans paralleles
			if((p1.getB()==0)&&(p2.getB()!=0)){
				Plane temp=p1; p1=p2; p2=temp;
			}
			// Chercher un point (x,y,z)
			double bdm=p2.getB()*p1.getD()-p1.getB()*p2.getD(); 
			double cbm=p2.getC()*p1.getB()-p1.getC()*p2.getB(); 
			double zconstante=bdm/cbm;
			double yconstante=-p1.getC()/p1.getB()*zconstante-(p1.getD()/p1.getB()); 
			Pos3D lepoint=new Pos3D(0,yconstante,zconstante); 
			Pos3D ladirection=new Pos3D(1,0,0);
			return new DroiteSpace3D(lepoint,ladirection); 
		}
		
		double ab=p1.getA()*p2.getB()-p1.getB()*p2.getA(); 
		if(ab==0) // vecteur directeur dans le plan xOy (composante nulle selon z)
		{
			double bdm=p2.getB()*p1.getD()-p1.getB()*p2.getD(); 
			double cbm=p2.getC()*p1.getB()-p1.getC()*p2.getB(); 
			double zconstante=bdm/cbm;
			double yconstante=-1/p1.getB()*(p1.getC()*zconstante-p1.getD());
			Pos3D lepoint=new Pos3D(0,yconstante,zconstante);
			Pos3D ladirection=new Pos3D(1,-p1.getA()/p1.getB(),0); 
			return new DroiteSpace3D(lepoint,ladirection); 
		}
		double adm=p2.getA()*p1.getD()-p1.getA()*p2.getD(); 
		double acm=p2.getA()*p1.getC()-p1.getA()*p2.getC();
		double yconstante=adm/ab; 
		double zconstante=0; 
		double xconstante=-p1.getB()/p1.getA()*(adm/ab)-p1.getD()/p1.getA(); 
		double zcoef=1; 
		double ycoef=acm/ab; 
		double xcoef=-p1.getB()/p1.getA()*(acm/ab)-p1.getC()/p1.getA(); 
		Pos3D lepoint=new Pos3D(xconstante,yconstante,zconstante);
		Pos3D ladirection=new Pos3D(xcoef,ycoef,zcoef); 
		return new DroiteSpace3D(lepoint,ladirection); 
	}
	
		/* Trouver l'intersection de la droite avec un des quatres cotes de l'image
		 * retourne null si parallele ou en dehors de l'image. 
		 * on ne s'occupe pas de la valeur de x (tout se passe dans le plan perpendiculaire
		 * a l'axe Ox)
		 */
		public Pos3D coteVsDroite(Pos3D A,Pos3D B,DroiteSpace3D D){
			double yv=D.getDirection().getY(); 
			double zv=D.getDirection().getZ(); 
			double ya=A.getY(); 
			double za=A.getZ(); 
			double yb=B.getY(); 
			double zb=B.getZ(); 
			double y0=D.getPoint().getY(); 
			double z0=D.getPoint().getZ(); 
			double det=-zv*(ya-yb)+yv*(za-zb); 
			if(Math.abs(det)<1e-6) return null; 
			double detAlpha=-zv*(y0-yb)+yv*(z0-zb);
			double alpha=detAlpha/det; 
			if ((alpha<0)||(alpha>1)) return null; 
			// on s'en fout de X, on va tracer une droite dans un plan
			double ny=alpha*ya+(1-alpha)*yb; 
			double nz=alpha*za+(1-alpha)*zb;
			return new Pos3D(0,ny,nz); 
		}
		
		public ArrayList<Pos3D> intersectCadre(double ymin,double ymax,double zmin,double zmax,DroiteSpace3D Q){
			ArrayList<Pos3D> resu=new ArrayList<Pos3D>(); 
			Pos3D A=new Pos3D(0,ymin,zmin); 
			Pos3D B=new Pos3D(0,ymin,zmax); 
			Pos3D C=new Pos3D(0,ymax,zmin); 
			Pos3D D=new Pos3D(0,ymax,zmax); 
			if(coteVsDroite(A,B,Q)!=null) resu.add(coteVsDroite(A,B,Q)); 
			if(coteVsDroite(A,C,Q)!=null) resu.add(coteVsDroite(A,C,Q)); 
			if(coteVsDroite(B,D,Q)!=null) resu.add(coteVsDroite(B,D,Q)); 
			if(coteVsDroite(C,D,Q)!=null) resu.add(coteVsDroite(C,D,Q)); 
			return resu; 
		}
	
	public void run(String arg) {
		
		// Les parametres de la cyclide
		double a=5;
		double mu=3; 
		double c=1; 
		
		// le centre du cercle
		double cx,cy,cz; 
		
		// le rayon du cercle
		double rayon; 
		
		// Les infos du plan
		double v1=0,v2=0,v3=0,pscal=0; 
		
		primeIP = new ColorProcessor(X,Y);
		
		
		double limit=1; 
		double lmoins=-1;
	
		setBounds(lmoins,limit,lmoins,limit); 
		setCoefs(); 
		
		
		
		Polygon background=new Polygon(new int[]{0,0,X,Y},new int[]{0,X,Y,0},4);
		primeIP.setColor(Color.WHITE);
		
		primeIP.fillPolygon(background);
		
		
		primeIP.setFont(new Font(null,Font.ITALIC,60));
		primeIP.setColor(Color.BLACK);
		primeIP.setLineWidth(2); 
		

		setProcessor("Verifying", primeIP);
		
		
		primeIP.setColor(Color.BLACK);
		primeIP.setLineWidth(5);
		Mesh leMesh=null; 
		try{
		leMesh=STLReader.readMesh("/tmp/cyclide.stl");
		
		BufferedReader in = new BufferedReader(new FileReader("/tmp/planesExtended.txt")); 
		
		double ymin=100; 
		double zmin=100; 
		double ymax=-100; 
		double zmax=-100; 
		//int nbCercles=4; 
		Ligne memoire[]=new Ligne[nbCercles]; 
		Pos3D lesCentres[]=new Pos3D[nbCercles]; 
		//Circle3D lesCercles[]=new Circle3D[nbCercles];
		Circle3D lesCerclesAvant[]=new Circle3D[nbCercles];
		for(int i=0;i<nbCercles;i++){
	    String ligne = in.readLine();
       
        Scanner rl=new Scanner(ligne); 
		cx=rl.nextDouble(); 
		cy=rl.nextDouble(); 
		cz=rl.nextDouble(); 
		
		Pos3D imCenter=new Pos3D(cx,-cz,cy); 
		
		
		rayon=rl.nextDouble();
		
		
		
		v1=rl.nextDouble(); 
		v2=rl.nextDouble(); 
		v3=rl.nextDouble(); 
		pscal=rl.nextDouble(); 
		
		
		
		Plane lePlan=new Plane(v1,-v3,v2,pscal); 
		System.out.append("Plan lu "); 
		System.out.println(lePlan); 
		
		Circle3D leCercle=new Circle3D(imCenter,rayon,lePlan); 
		lesCerclesAvant[i]=leCercle; 
		lesCercles[i]=leCercle.rotate(rotations(lePlan));
		System.out.println("Cercle :"+lesCercles[i]);
		
		Ligne cochonou=tranche(leMesh,lePlan);
		cochonou=cochonou.rotate(rotations(lePlan)); 
		memoire[i]=cochonou; 
		lesCentres[i]=imCenter.rotate(rotations(lePlan)); 
		
		System.out.println(cochonou.info()); 
		if(cochonou.getMinY()<ymin) ymin=cochonou.getMinY(); 
		if(cochonou.getMinZ()<zmin) zmin=cochonou.getMinZ(); 
		if(cochonou.getMaxY()>ymax) ymax=cochonou.getMaxY(); 
		if(cochonou.getMaxZ()>zmax) zmax=cochonou.getMaxZ(); 
		
		

		//print(cochonou,20,0.1,"wtest",i); 
		}// for
		
		intersections=new DroiteSpace3D[nbCercles][nbCercles]; 
		for(int i=0;i<nbCercles;i++){
			Pos3D pointCourant=pointCommun.rotate(rotations(lesCerclesAvant[i].getPlan()));
			for(int j=0;j<nbCercles;j++)
			{
				intersections[i][j]=DroiteSpace3D.inter(lesCerclesAvant[i].getPlan(), lesCerclesAvant[j].getPlan()); 
				System.out.println("\n ***********************\n"); 
				System.out.println("Plan 1 :"+lesCerclesAvant[i].getPlan()); 
				System.out.println("Plan 2 :"+lesCerclesAvant[j].getPlan()); 
				System.out.println("Droite obtenue par intersection : \n"+intersections[i][j]+"\n"+intersections[i][j].belongs(pointCommun));
				System.out.println(pointCommun);
			}
		}
		
		System.out.println("Y min : "+ymin+" Y max : "+ymax); 
		System.out.println("Z min : "+zmin+" Z max : "+zmax); 
		double spreadMax; 
		if((ymax-ymin)>(zmax-zmin)) spreadMax=ymax-ymin; else spreadMax=zmax-zmin; 
		for(int i=0;i<nbCercles;i++){
			if(i<nbCercles/2)
			print(memoire[i],lesCerclesAvant,spreadMax,0.1,"AAtest",i,true); 
			else{
				print(memoire[i],lesCerclesAvant,spreadMax,0.1,"BBtest",i,false); 
			}
		}
		
		
		}//try
		catch(Exception e){System.out.println(e); System.exit(0); }
		
		
		System.out.println("fini"); 
		
		
	
}
	
}
	